//
//
//      +-----------+-----------+
//      |           |           |
//      |         (i,j)         |
//      >     *     >     *     >
//      |   (i,j)   | (i+1,j)   |
//      |           |           |
//      +-----------+-----------+
//                       
//                  |           |           |           |
//                --^-----------^-----------^-----------^-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |  (i,j+1)  |     :     |
//                  o     >     o    u_N    o     >     o   
//                  |     :     |     :     |     :     |
//                  |     :     |     :     |     :     |
//                --^---------- 3 -- v_n -- 4 ----------^--  4 = v(i+1, j  )
//                  |     :     |     :     |     :     |    3 = v(i  , j  )
//                  |     :     |     :     |     :     |    2 = v(i+1, j-1)
//                  o    u_W   u_w   u_P   u_e   u_E    o    1 = v(i  , j-1)
//                  |  (i-1,j)  |   (i,j)   |  (i+1,j)  |
//                  |     :     |     :     |     :     |
//                --^---------- 1 -- v_s -- 2 ----------^-- 
//                  |     :     |     :     |     :     |               
//                  |     :     |     :     |     :     |
//                  o     >     o    u_S    o     >     o   
//                  |     :     |  (i,j-1)  |     :     |
//                  |     :     |     :     |     :     |
//                --^-----------^-----------^-----------^--
//                  |           |           |           | 
//                   
//   u_w = ( u(i-1,j) + u(i,j) ) / 2     u_e = ( u(i+1,j) + u(i,j) ) / 2
//   v_n = ( v(i,j) + v(i+1,j) ) / 2     v_s = ( v(i,j-1) + v(i+1,j-1) ) / 2
//              3        4                            1          2          
// 
//---------------------  3D  ---------------------
//

namespace Tuna {

template<class Tprec, int Dim>
inline
bool CDS_XLES<Tprec, Dim>::calcCoefficients(const ScalarField &nut) {
    Tprec dyz = dy * dz, dxz = dx * dz, dxy = dx * dy;
    Tprec dyz_dx = dyz / dx, dxz_dy = dxz / dy, dxy_dz = dxy / dz;
    Tprec ce, cw, cn, cs, cf, cb;
    Tprec nutinter;
    Tprec dxyz_dt = dx * dy * dz / dt;

    for (int i =  bi; i <= ei; ++i)
	for (int j = bj; j <= ej; ++j)
	    for (int k = bk; k <= ek; ++k)
	    {
		ce = ( u(i+1, j, k) + u(i,j,k) ) * 0.5 * dyz;
		cw = ( u(i-1, j, k) + u(i,j,k) ) * 0.5 * dyz;
		cn = ( v(i,j,k) + v(i+1,j,k) ) * 0.5 * dxz;
		cs = ( v(i,j-1,k) + v(i+1,j-1,k) ) * 0.5 * dxz;
		cf = ( w(i,j,k) + w(i+1,j,k) ) * 0.5 * dxy;
		cb = ( w(i,j,k-1) + w(i+1,j,k-1) ) * 0.5 * dxy;
//
// nut is calculated on center of volumes, therefore, nut
// must be staggered in x direction:
		nutinter = 0.5 * ( nut(i,j,k) + nut(i+1,j,k) );
		
		aE (i,j,k) = (Gamma + 2 * nutinter) * dyz_dx - ce * 0.5;
		aW (i,j,k) = (Gamma + 2 * nutinter) * dyz_dx + cw * 0.5;
		aN (i,j,k) = (Gamma + nutinter) * dxz_dy - cn * 0.5;
		aS (i,j,k) = (Gamma + nutinter) * dxz_dy + cs * 0.5;
		aF (i,j,k) = (Gamma + nutinter) * dxy_dz - cf * 0.5;
		aB (i,j,k) = (Gamma + nutinter) * dxy_dz + cb * 0.5;
		aP (i,j,k) = aE (i,j,k) + aW (i,j,k) +	
		             aN (i,j,k) + aS (i,j,k) +	    
		             aF (i,j,k) + aB (i,j,k) + 
		             dxyz_dt;
//		aP (i,j,k) /= alpha;  // under-relaxation
//		+ (ce - cw) + (cn - cs) + (cf - cb);	    
// Term (ce - cw) is part of discretizated continuity equation, and
// must be equal to zero when that equation is valid, so I can avoid
// this term for efficiency.

		sp (i,j,k) =  u(i,j,k) * dxyz_dt -
		    ( p(i+1,j,k) - p(i,j,k) ) * dyz +
		    nutinter * ( (v(i+1,j,k) - v(i,j,k) - 
				  v(i+1,j-1,k) + v(i,j-1,k)) * dz +
				 (w(i+1,j,k) - w(i,j,k) - 
				  w(i+1,j,k-1) + w(i,j,k-1)) * dy );
		  
//		    u(i,j,k) * (1-alpha) * aP(i,j,k)/alpha;// under-relaxation
	    }

    calc_du_3D();
    applyBoundaryConditions3D();

    return 1;
}

} // Tuna namespace

















